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We develop the multiscale technique to describe excitations of a Bose- Einstein condensate (BEC) 
whose characteristic scales are comparable with the healing length, thus going beyond the conven- 
tional hydrodynamical approximation. As an application of the theory we derive approximate 
explicit vortex and other solutions. The dynamical stability of the vortex is discussed on the basis 
of the mathematical framework developed here, the result being that its stability is granted at least 
up to times of the order of seconds, which is the condensate lifetime. Our analytical results are 
confirmed by the numerical simulations. 

PACS number(s): 03.75.-b, 03.75.Fi, 67.57.Fg 



d . I. INTRODUCTION 

a 

' Since the remarkable experimental realization |jj of the Bose-Einstein condensation (BEC) there has been an 
explosion of the experimental and theoretical activity devoted to the physics of dilute ultracold bosonic gases Q . 

From the theoretical point of view, many of current condensate systems are cold enough that the mean field theories 
can be applied to describe the ground state of the system S. The most classical model arising in this context for 
the order parameter of the system is the the Gross-Pitaevskii (GP) equation. In fact, many relevant features of 
' the condensate such as free expansion B, low energy excitations p], collapse HQ], dynamics of multicomponent 
. systems ||], vortex generation and dynamics ]9|,[l0|, and many others, can be understood within the framework of this 
CO ' "simplest" mathematical model. 

Mathematically, the GP equation, to be written explicitly below, is an equation of Nonlinear Schrodinger type. 
These equations have been studied extensively both in the physical and mathematical literatures []l2|,[l3| since they 
provide an universal model for the study of the dynamics of envelope waves. Some fields where these equations 
arise are plasma physics [OL fiber and integrated optics [ fl5f , wave propagation in Kerr media |l6f| , water waves 



OO 
O 

| p7[ , biomolecule dynamics |18[ and fundamentals of quantum mechanics 19 . One of the distinctive features of the 
equation as it appears in BEC problems is the presence of an external potential -the trap-, which affects essentialy 
d ' the spectrum of the "linear part" , i.e. of the underlying Schrodinger equation. 

The mathematical analysis of the GP equation in the situations common in BEC is very complicated and this is why 
most theoretical studies deal with either numerical simulations or some kind of approximations such as the Thomas- 
Fermi, the hydrodynamic approach or simply reduce the number of degrees of freedom by using the variational 



approach based on trial functions. To be precise let us write explicitly the equation for the order parameter, whose 
form is 

ih^ = -^V 2 * + V ext (r)* + g\*\ 2 *, (1.1) 
where the external potential can be well approximated by 



V ext = ^m(LJ 2 x x 2 + Ljy+u; 2 z z 2 ), (1.2) 
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and g = AirTi 2 a s /m. 

The simplest approximation that can be done to obtain information on the ground state is the so called Thomas- 
Fermi approximation, which corresponds to neglecting the kinetic energy term under the assumption that it is much 
less than the self-interaction energy. In the cases when the trap is smooth enough this approximation is not bad 
although it is not consistent, since the assumption of zero kinetic energy provides a profile that has an infinite value 
of that quantity. 

The following level of approximation, the hydrodynamic approach, can be explained by introducing new variables 
n(r, t) and v(r, t) through the relations 

#(r,t) = V^Oe 15 ^, (1.3a) 
v(r,t) = — WS(r,t). (1.3b) 

777 

Then, one arrives to the system of equations 

d t n + V • (vn) = 0, (1.4) 



/ 7j TfllP' \ 

md t v + X7[^V ext +gn- 1 ^^V 2 V^+—)=0. (1.5) 

The conventional hydrodynamic approximation corresponds to neglecting the term proportional to ^7 2 y/n under the 
assumption that the density is a slowly varying function p3| . For a spherically symmetric potential (uj x — uj v = oj z = 
ui) this can be done in the limit Na s /ao 3> 1, where a s > is the scattering length and ag = y/h/mu) is the harmonic 
oscillator length, i.e. the trap size. 

To the best of authors knowledge the mathematical limits of validity of the above mentioned approximations have 
not been established (which would imply study of the next orders of the expansion or making appropriate rigorous 
estimates), so far. In particular, using the estimate n ~ 10 12 cm -3 one finds that for a s ~ 10 nm the term with 
S7 2 y/n/ \fn is much less than gn for excitations whose characteristic length A is of order of 10 -3 cm. Such a A is 
comparable with the characteristic size ao of the trapping potential. Thus the question about analytical description 
of the excitations which vary on scales much less than ao appears. It is important to notice that in those cases the 
terms ^ V /^r"' arL< ^ 9 n become comparable with each other and have to be treated in equivalent manner. Using the 
wave terminology this situation can be characterized as balance of dispersive and nonlinear effects, which should lead 
to stable spatially localized excitations. 

It is our intention in this paper to describe an analytical tool, the multiscale approximation, which is able to provide 
approximate solutions of the problems under consideration in explicit analytical form. In less rigorous terms this 
technique has been applied in the context of BEC in Refs. j7pC|] to problems with negative scattering length. The 
application to positive scattering length presents some new characteristics as we will see later (in particular it implies 
different scaling order). Specifically, we will consider a BEC in a parabolic cilindrically or spherically symmetric 
confining potential and obtain a set of background solutions in the form of black holes or vortices whose characteristic 
spatial scale is much larger than the range of inter-atomic interactions but smaller than the trap size. Working in this 
way we will be able not only to obtain an approximate analytical shape for the vortex solutions but also to provide 
results on their dynamical stability, a point which is still of great interest after the many studies devoted to it |^l|,^2| . 
The outline of the paper is as follows. First, in Sec. |l| we pose the problem in precise terms. As a second step, we 



apply the technique of multiscale expansion to our particular problem in Sec. Ill 



set of solutions and discuss their stability. Finally, in Sec. |y|, we summarize our conclusions. 



In Sec. IV we obtain a particular 



II. STATEMENT OF THE PROBLEM 

We start with the three-dimensional GP equation in the form 

ifi^ = -^-V 2 ^ + V ext {r)^ + J A (r)*. (2.1) 
We will preserve the integral form of the interaction term [c.f. Eq. ( |1.1| )] 

J X (r)=g I T/ A (r,r')|*(r',i)| 2 dr', (2.2) 
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since there are situations where its role could be important [£3| and it does not complicate the formalism. Here 
V\(r,r') describes two-body interactions and A is the effective interaction length (mathematically, the size of the 
kernel) . 

The wave function is normalized to the total number of particles N 

J |*(r)| 2 dr = N. (2.3) 

It is reasonable to assume that the effective interaction length A parametrizes the kernel so that 

limV A (r,r') =S(r-r'). (2.4) 

The last limit corr esponds to the usual case of local interactions and gives rise to the |\E'| 2 nonlinear term in the GP 
equation [i.e. (2.1) transforms into fll-l| )]- Bearing in mind this last fact we can narrow the class of the kernels under 
consideration by requiring 

Vi(tJ) = ±u(^^), (2.5) 

where 

J U(\i -?])<& = 1. (2.6) 
and f = r/A. As it is clear from the definition of A, it must be of order of the scattering length: 

- = 0(1) (2.7) 

a s 

In what follows, however we keep A as a parameter. 

We will concentrate on a BEC in a trap with either cilindrical (with z-axis being the symmetry axis) or spherical 
symmetry, i. e. oj a = u, where a = x, y in the former case and a = x, y, z in the last case. We define the characteristic 
scales of the condensate in the trasversal, ao = yjh/(muj), and longitudinal, a z — ^Jh/(muj z ) 1 directions. In the case 
of spherical symmetry these scales coincide, a z = ao. 

In order to provide a formal description of the system it is convenient to scale out the GP equation. To this end 
we introduce the dimensionless time and wavefuntion through 

t, (2.8a) 



2mA 2 



#Sr)=2Ay^*(r,t). (2.8b) 
Then the renormalized GP equation takes the form 

M = -Aij + V ext {t)^ + j(i)^ (2.9) 

OT 

where 

J(f)= /[/(f-f')|^(f',r)| 2 df', (2.10) 



V ext (r) 



(Sl 2 f ^ + S! 2 .? 2 ) cilindrical symmetry 
fl 2 r 2 spherical symmetry 



(2.11) 



r_L = (x,y), O = A 2 /ag and O z = A 2 /a 2 . 

Having scaled out the kernel which describes the nonlocal interactions we have to impose the condition that the 
characteristic scale of the kernel variation is of order one, i. e. 



1 dU(f) 



O(l) 



U{f) dr 

where f = |f | and f must be substituted by fx in the case of the cilindrical symmetry. 
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III. MULTISCALE ANALYSIS. 



The multiscale technique is a mathematical technique for the analysis of problems where there are different spatial 
scales involved in the solution. The problem at hand possesses three important scales: the scattering length a s , 
the trap size ao and the healing length £ = (87ma s ) -1 / 2 . We will concentrate on the analysis of excitations of the 
superfluid corresponding to the scale of the healing length, which occur on top of a background with spatial scale 
of the order of the trap size. In other words our theory will be valid for cases when the above three parameters are 
related as follows 



a s <£<a . (3.1) 
For the first step we recall the properties of the linear spectral problem 

L d \n,l,m) = £ n ,i,m\n,l,m), (3.2) 
where = Cd + M.d- The operator Cd describes the radial dependence 



P i- 1 dp P0 d P0 



C * = - — —jt— + n Po, (3-3) 

and 



Q2 X Q2 

M * = 7^2 + 72 TTZa, ( 3 - 4a ) 



dz 2 Pq d<t> 2 ' 
_L d_ . d_ 1_JP 

p 2 o sm0d9 Sln 86 + p 2 d(j) 2 



-^3 = 72— 7T^ sin ^ + T2^l2- (3-4b) 



The indexes n and m refer to radial and magnetic quantum numbers while I is the azimutal (in the spherical case) or 
longitudinal (in the cilindrical case) quantum number (see Appendix) . Hereafter d = 2 and d = 3 are associated with 
cilindrical and spherical symmetries. 

For the next consideration it is convenient to single out the radial part of the eigenfunctions and represent 

I 7 \ J £n,m (po)O,m(z,0) cilindrical case ,„ ^ 

' ' \€n,l(Po)(l,m(0,(t>) spherical case 1 ' ; 

To shorten notations in what follows we drop the second subindex of £ in the cases at hand: = £„, keeping in 
mind that n is a radial (rather than total) quantum number. 

The operator Cd is considered in the space of functions satisfying the conditions 

|£ n (0)| < oo, lim po£n(po) = 0. 
The inner product in this space is defined by 

poo 

(£n,fn')= / Up)^{p) P d - 1 dp. (3.6) 



Then the operator C is Hermitiam. 

The problem at hand is directly related to the well known eigenvalue problem for the linear oscillator in two and 
three dimensional parabolic potential p5| . Its spectrum is discrete 

_ J Ail (n + \ + y) + 2f2 z (l + |) cilindical case , , 

^-140 (»+| + {) " spherical case 

(both cases naturally coincide after substitution il z i— > VL and I + m <— > I in the first formula). The orthormalized 
eigenfunctions of Cd are given in the Appendix. 

In accordance with the multiscale method, which is a theory for weakly nonlinear problems, the solution is searched 
in the form of a weakly modulated linear mode. This is a main feature of the multiscale method, which means that 
the theory is applicable to small condensates. However, in many problems the results coming out from the multiscale 
expansions remain valid even for t he l arge nonlinearity limit. Respectively we introduce a formal small parameter 
e -C 1 and look for the solution of (2.9) in the form 
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tp = eipi + e ip 2 + 



(3. 



Some peculiarities of the expansion introduced are to be discussed here. First, we use the expansion of only radial 
part of the solution. This means that the consideration will be restricted to symmetric solutions (the generalization 
is straightforward but cumbersome). Second, in what follows we will concentrate on modulation of the ground state 
which turns out to be flat compared with the excitation itself. This however is not in contradiction with the assumption 
about smooth (i.e. dependent on pi) modulation of the underlying linear mode, depending on the rapid variable po, 
by weak nonlinearity. This is justified by the following reasons: (i) The background solution must possess the "right" 
assymptotic behavior at p — > oo, which can be obtained only by taking into account the term fl 2 p 2 at the "linear 
level" ; (ii) The expansion must take into account the contribution of all linear eigenmodes, rather than only the 
leading one [see ( 3 . 13j ) and discussion below] ; (iii) The linear problem with the parabolic trap is solved exactly. 



The functions sJP~ depend on the set of spatial coordinates pj = e J p and times Tj = e>T (j = 0, 1,2, ...): ^ = 

, px, p2, ...;to,ti, T2, ...). Hereafter p = \r±\ in the cilindrical case and p — |f| in the spherical case. The variables 
pj and Tj are regarded as independent: thus 



1 d 



P' 



d-1 



dp 



p 



d_ 

dr 

dp 



d_ 

dr Q 
1 



d 



d 



d d-i 9 
~dp~o P ° dio + 



Po 



1 d 



d p d pi 



Po dpi 



0(e 2 



(3.9a) 
(3.9b) 



Before going into details of the multiscale expansion we have to clarify the physical meaning of the parameter e. 
We are interested in excitations against a background. Then in terms of the dimcnsionless wave function tpi the 

where V 



normalization condition (2.3) leads to the estimate N 



<\ 2 



is the volume of the condensate 



which is limited by the confining potential. These estimates imply the normalization condition 



d n 3-d 



\M^r)\ 2 dr 



(3.10) 



(notice that in the case of cilindrical symmetry, taking A = a s one obtains that the new wave function is normalized 
to one) . In other words the effective small parameter of the multiscale expansion can be identified as 

(3.11) 



where n — N/V is the density of the condensate and (2/7) is used. Taking into account that typically na^ — 
O(10~ 4 -j- 10~ 6 ) one concludes that e (~ 1CP 2 -j- 1CP 3 ) can indeed be used for the multiscale expansion. We notice 
that the approximation of the local kernel, A — > is evidently a particular limit of the theory. Formally, however, it 
is more convenient to keep e as a small parameter without specification and substitute it by one in the final formulas. 

We are interested in the evolution of an initially unperturbed "linear" state [say, (no, Iq, mo)th one]. Gathering the 
terms of the first order with respect to e expansion we conclude that ip\ can be represented in the form (£ = £ n0y i Q . m o) 



d 1} - Mpi 



;n, . 



- i£ " T "\n ,l ,m ). 



(3.12) 



Hereafter the notation (p 1? t±, ...) is used in order to indicate the dependence on all spatial and temporal variables, 
which are slower than those written explicitely [i.e. for example in ( [3.12| ) the amplitude depends on t\,T2, ■■■■]■ 
The solution of the equation appearing in the second order of e is represented in the form of the expansion 



V^ 2 = J2 B n'(pi,-;ri,...)e- i£oTo \n,l,m) 



(3.13) 



Hereafter a prime means that the sum is computed over all arguments such that n ^ uq, I ^ Iq, and m ^ tuq. The 
respective equation of the second order takes the form 



y7(£p _ £n,l,m)B n ,l,m\n, I, m) + 



0A n 



d d-l\ dA n 



\no,lo,rn ) = 



(3.14) 



dn V dPo Pa J dpi 

Applying (no, Iq, mo\ to this equation (see Appendix) one arrives at the relation dA/dri = 0. Here we use the property 
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Kn,#n) = * / £nM [2^- + ^— ^ ) £„(p )Po ^0 = 

o V "Pa p a 



with 



P = *(^ + —) (3-15) 
V d Po Pa ) 

being the operator of the radial component of the (dimensionless) linear momentum. Then one concludes that A does 
not depend on t\, A = A[px, r 2 , ...). 

The coefficient B nn j n , m n is equal to zero identically [otherwise it could be made zero by simple renormalization of 



initial conditions for (2.E)]. In order to find other coefficients B n ,;,m we compute the inner product with (rio,la,Tno\. 



Introducing matrix elements 

Pn'n = {(n>P€n>) (3.16) 

corresponding to the transition between the states n' and n, originated by the radial component of the linear momen- 
tum considered as a perturbation, we obtain the coefficients of the second order term 

B n ,i,m(pi, •■■; T "2, ■•■) = i c Pn ° n (3.17) 

~ £-n,l,m OPl 

Finally, collecting all the terms of the third order of e and using the explicit form of the second order addendum 
found above we arrive at the dynamical equation for the amplitude A 

dA d 2 A 

t d^ +V °^p2-^ono\A\ 2 A = Q (3.18) 



Here 



is the effective dispersion and 



v Q = v no , , mo = i - Y! p lPn ° n? f ( 3 - 19 ) 



ii .1 . m 



e n {p)^(p')U{\i-~v'\)didv' (3.20) 

is the effective nonlinearity. As a matter of fact the effective nonlinearity can be computed under the assumption 
that the characteristic scale of t he g ro und state is much lar ger than the width of the kernel U(\r |) of the interaction 
potential. Taking into account (2.7), ( |3.l| ), and (jAll), (A3) one concludes tha t thi s supposition is valid with great 



accuracy for several lowest eigenstates £ ra (p). Then £/(|r|) in the integrand of ( p. 20] ) can be substituted by the delta 
function and thus \ nn — 

In the explicit form of T>o one can recognize the second order addendum to the energy of an effective particle 
induced by the perturbation operator p. Thus the effective dispersion has a contribution of the probability of the 
energy transfer among the levels in the confining potential. Notice that the mentioned transition occur among levels 
with the sa me orb ital and azimutal quantum numbers (in the case at hand they are equal to zero). 



Equation (3. IS) is considered on the semiline (p\ > 0) subject to the boundary condition lmi pl ^ oc A(p\, T2, ...) — 
Ae luJT2 where the constant A can be interpreted as an amplitude of the fundamental state of the condensate: after all 
it must be found from the normalization condition. [Notice that other boundary conditions, of less practical interest, 
are also compatible with the problem at hand]. 

Let us now return to the discussion of the main scales of the problem. In the dimensionless variables the spatial 
scale of the amplitude A(pi, ...;T2, ...) can be estimated as 1/e. On the other hand the problem possess another 
important dimensionless scale 1 / which determines the wave function localization due to the confining potential. An 
interesting and physically relevant situation appears when e ^> y/Q, what in physical variables means 87r7Va s /ao ^> 1 
psf , a condition which is satisfied in current experimental setups. The last requirement can be understood equivalently 
as smallness of the healing length £ compared with the linear oscillator length £ <C ao- I n this situation the excitation 
width is much smaller than the volume "available" for the condensate and in the case of positive scattering length 
one can obtain localized excitations against the background. This is the regime of validity of our theory as has been 
put in more mathematical terms in the first part of this section. 
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IV. APPLICATION: VORTEX SHAPE AND STABILITY 



One of the simplest ex citati ons of the ground state condensate one can construct have the form of static " dark" 
soliton solutions of Eq. (3.18), i. e. "black holes" or vortices (when a local phase is added). Joining the different 
contributions from the multiscale method we obtain 



T 1 a s 



2 1 IJT\ ht 

— + — + 8irJ— ) — 



exp 



2a 2 , 



2a? 



tanh 



im<fi 



0(nal 



(4.1) 



This solution approximates the shape of a vortex line with axis of rotation z in trap of transverse size Oo- 

This is not the only solution one can construct, as an example we present also a solution corresponding to a "black 
hole" in three dimensions 



1 



TT 3 / 4 f 



exp 



. 2 r/s 



ht 

2m 



exp 



2al 



tanh 




0(na 3 s ) 



(4.2) 



where r = (x, y, z). 

In the above formulas we have neglected the difference between A and a s [see (2.7)]. They however still contain the 
quantity Vq to be estimated. To th is end we employ (Al) and ( |A3| ) and compute (3.19) p^ |. The result is Vq w 0.586 
and I?o ~ 0.728 for (4.1) and (4.2) respectively. Thus in the physical units the effective width of the black hole can 
be estimated to be of order of the healing lengths which is the expected result for the vortex size. 

It must be said that the result naturally does not depend on 17, a fact which is explained by taking into account 
that in the parabolic trap the energy levels are equidistant. Another point to be m entioned is that as a matter of fact 
we have computed the terms of order of 0(na 3 ). They are given by (3.13), (3.17) and are not represented explicitly 
because of cumbersome form. In other words the theory provides one with the BEC wavefunction computed with 
accuracy 0((na 3 ) 3 / 2 ) to problems whose number of particles is not too large. 

One of the advantages of the multiscale method is that it guarantees the dynamical stability of the solution up to 
a certain time scale. In order to explain this last point we notice that the method is based on eliminating the secular 
terms in the third order with respect to the small parameter. This in particular means that the characteristic time of 
evolution of instabilities of the solutions (if any) is given by e 3 r. The main implication is that the instabilities a gains t 
the background solution (which is the vortex in the case at hand) cannot develop up to r ~ e~ 3 . Recalling now (3.11) 
and returning to the physical variables we obtain that for the typical values of e ~ 10~ 2 any instability could develop 
only for times larger than 10 seconds, which is of the order of the life-time of the condensate. In fact, this prediction 
is consistent with the previous numerical observation |2^| that the vortex is a stable object. 

This result is relevant since a mathematical proof of the expected dynamical stability of a vortex in a trap (i.e. 
when dissipation is zero) is not available yet. However, the numerical simulations (which are also valid only for limited 
times) point that a vortex should be a dynamically stable object. The analysis of the linearized operator around the 
vortex done in Ref. |p2f was not completely conclusive since there were some zero eigenvalues which could lead to 
instabilities mediated by high order terms in the linearized evolution operator. In this sense the multiscale technique 
adds a new argument in favor of long time stability of the vortex in the framework of the mean field modelization of 
the phenomenon. 

We remark that the previous discussion applies to dynamical stability, i. e. stability of the solution under pertur- 
bations on the initial data. It was discussed on the previous works on vortex stability that the system is energetically 
unstable |2l],^2|], while the question of its Lyapunov stability remains open. Our analysis concerns only the subject 
of dynamical stability and confirm the numerical observations by many authors that vortices in a trap are very long 
lived and could even be completely stable p2|. 



V. CONCLUSION 



In the present paper we have developed the multiscale approach for obtaining explicit (approximate) solutions to 
the GP equation for condensates with sufficiently small number of particles which are characterized by fast spatial 
scales of the order of the healing length, and thus cannot be described within the framework of the conventional 
hydrodynamical approach. One of the advantages of the method is that it guarantees the stability of the solution 
as discussed in Sec. |iy] which allows us to put lower bounds on the time a vortex survives under possible dynamical 
instabilities described within the framework of the GP equation. 

Another outcome of the theory developed above is that the solutions of the healing length size cannot be excited 
only as a background state. Some accompanion modes, corresponding to higher levels of the confining potential are 
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excite d simultaneously. All these contributions are included in the constants T>q appearing in the reduced equations 



(M3 



One more advantage of the multiscale method is that, by reducing the problem to solving the effective one dimen- 
sional nonlinear Schrodinger equation, one gets a powerful tool for obtaining other more complicated solutions and 
even interactions among them, a point which will be the subject of further study. 

The present paper provides a new analytical tool which we hope will be of interest for theoretical progresses in the 
analysis of BEC problems. 
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APPENDIX A: EIGENFUNCTIONS 

For the sake of convenience here we present the orthonormalized eigenfunctions. In the case of potential with 
cilindrical symmetry they are 



£n,m(p) 



'2fi 



(n + m) 



(Al) 



1/4 



7T 3 / 4 V2 1 + T l! 



Here H n (-) is the Hermite polinomial and L„ (•) is the Laguerre polinomial. When spherical symmetry is imposed 
the eigenfunctions are given by 



„im4, e -n,t>/2 H ( ^ 



(A2) 



£n,l(p) = 2 



(2(n + + 



r(n+l) 



1/2 / 



(A3) 



Y, 



(-1) 



m+\m\-l 



21 + 1(1- \m\ 



4tt (Z + |m|)l 



1/2 



^'"'(cos^e 1 ™*. 



The inner product in this space is defined by 

/OO /•2-7T 
dz / i 
-oo JO 

in the case of cilindrical symmetry and by 



f7T /■27T 

(n, l,m\n',l',m') = \ d0 dcf) / dppZ, n , m {p)Z,n' ■,m'{p)Cl,m{OA)Cv 

7T JO 



(A4) 



(A5) 



(A6) 



in the case of spherical simmetry 
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